import sympy as sp

# 定义符号变量
x, y = sp.symbols('x y')

# 定义标量函数 φ = (4-x²)²(4-y²)²
phi = (4 - x**2)**2 * (4 - y ** 2) ** 2

# 计算旋度 curl(φ) 在二维中定义为 (dφ/dy, -dφ/dx)
u1 = sp.diff(phi, y)  # 旋度的x分量
u2 = -sp.diff(phi, x)  # 旋度的y分量

# 定义向量场u
u = (u1, u2)
print("向量场u(x, y)的表达式:")
print(f"u1 = {sp.simplify(u1)}")
print(f"u2 = {sp.simplify(u2)}")

# 计算拉普拉斯算子 Δu = (Δu1, Δu2)，其中Δui = d²ui/dx² + d²ui/dy²
delta_u1 = sp.diff(u1, x, 2) + sp.diff(u1, y, 2)
delta_u2 = sp.diff(u2, x, 2) + sp.diff(u2, y, 2)

# 简化结果
delta_u1_simplified = sp.simplify(delta_u1)
delta_u2_simplified = sp.simplify(delta_u2)

print("\n拉普拉斯算子Δu的表达式:")
print(f"Δu1 = {delta_u1_simplified}")
print(f"Δu2 = {delta_u2_simplified}")

# 进一步因式分解以获得更简洁的表达式
delta_u1_factored = sp.factor(delta_u1_simplified)
delta_u2_factored = sp.factor(delta_u2_simplified)

print("\n因式分解后的拉普拉斯算子Δu:")
print(f"Δu1 = {delta_u1_factored}")
print(f"Δu2 = {delta_u2_factored}")

### 补充计算梯度 grad(u) ###
# 计算u1对x和y的偏导数
du1_dx = sp.diff(u1, x)
du1_dy = sp.diff(u1, y)

# 计算u2对x和y的偏导数
du2_dx = sp.diff(u2, x)
du2_dy = sp.diff(u2, y)

# 构建梯度矩阵 (2x2雅可比矩阵)
grad_u = sp.Matrix([
    [du1_dx, du1_dy],
    [du2_dx, du2_dy]
])

# 简化梯度矩阵
grad_u_simplified = sp.simplify(grad_u)

# 对梯度矩阵的每个元素进行因式分解
grad_u_factored = sp.Matrix([
    [sp.factor(grad_u_simplified[0, 0]), sp.factor(grad_u_simplified[0, 1])],
    [sp.factor(grad_u_simplified[1, 0]), sp.factor(grad_u_simplified[1, 1])]
])

print("\n向量场u的梯度 grad(u) 矩阵:")
print("梯度矩阵（未简化）:")
sp.pprint(grad_u)

print("\n简化后的梯度矩阵:")
sp.pprint(grad_u_simplified)

print("\n因式分解后的梯度矩阵:")
sp.pprint(grad_u_factored)
